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Abstract 

• A coupled quantum-classical model describing the transport of electrons confined 

in nanoscale semiconductor devices is considered. Using the subband decomposition 
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approach allows to separate the transport directions from the confinement direction. 
The motion of the gas in the transport direction is assumed to be classical. Then a 
hierarchy of adiabatic quantum-classical model is obtained, leading to subband SHE 
" and energy-transport models, with explicit expression of the diffusion coefficients. The 

energy-transport-Schrodinger-Poisson model is then used for the numerical simulation 
of the transport of the electron gas in an ultra-scaled Double-Gate-MOSFET. 
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tor, Gummel iterations, mixed finite elements. 
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1 Introduction 



In nanoscale semiconductor devices, electrons might be extremely confined in one or several 
directions refered to as the confining direction. This leads to a partial quantization of the 
energy which can be modelled by the subband decomposition method [32]. This subband 
decomposition approach allows to separate the confinement direction from the transport 
direction. Thanks to the reduction of the dimension of the transport problem, the compu- 
tational gain is significant. In the confined direction electrons behave like waves; the system 
is at thermodynamical equilibrium and it is described by a statistical mixture of eigenstates 
of a Schrodinger-Poisson system. In the transport direction the transport can be of classical 
[8j [7] or quantum nature |32| . 

Here, we are interested in deriving adiabatic quantum-classical models accounting for 
thermal effects, aiming at accurate and efficient numerical simulation of confined devices. 
In [9] several spherical harmonic expension (SHE) models incorporating quantum effects 
are proposed. However, with their strategy the obtained models have a complicated non- 
local structure which is not suitable for numerical purposes. Quantum energy-transport 
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and quantum drift-diffusion models have been derived in [IS] using the strategy of quantum 
moments, as well as in (25]. These models involve a quantum chemical potential that depends 
on the density in a non-local way. In this work we follow the complementary strategy 
proposed in [7] where the subband model is derived first, then a diffusive approximation 
of the adiabatic Boltzmann equation is performed to obain coupled quantum-fluid models 
(spherical harmonic expansion and energy-transport). The subband energy-transport (ET) 
model in [7J, directly derived from the Boltzmann equations (as in e.g. j6] for the classical 
case), is, however, not immediately suited for numerical simulations, since the diffusion 
coefficients are not given in explicit form and, moreover, the energy relaxation term is 
not obtained. Therefore, we propose in this work a suitable description of the dominant 
collision mechanisms which allows to extend the formal derivation of the spherical harmonic 
expansion (SHE) model given in [7J. Then a new energy-transport (ET) model is formally 
derived as diffusive limit from this SHE model [3]. 

Numerical discretization of classical ET equations has already been studied in many 
papers : by using mixed finite elements schemes e.g. in [T5| [T9| [20| |22| [27J [29], ENO 
schemes in [23J, finite difference methods [T5| [M] and finite volume schemes in [13J. In 
[151122], the authors propose a drift-diffusion reformulation which allows to use an accurate 
Sharfetter-Gummel scheme with exponential fitting [11] and, moreover, to decouple the ET 
model. However, in this quantum case, the involved form of the diffusion coefficients does 
not allow for a decoupled drift-diffusion reformulation. Then we will use a more traditional 
approach with mixed finite elements which can be directly applied since the obtained ET 
model turns out to be in symmetric form. Then, a Gummel type algorithm is used as outer 
iterations of the solution of the coupled energy-transport-Schrodinger-Poisson model, and 
the (non-linear) ET discrete system is solved by means of a Newton scheme. Moreover, 
passing to the limit in the energy relaxation term, a subband drift-diffusion equation is 
recovered in the form of [5] with a more accurate description of the diffusion coefficients 
taking into account the collisional mechanisms (see also |31j for numerical simulations). 
This work is then an extension of [31] to a more general collisional framework. 

The outline of the paper is the following. In Section [2] we set the assumptions on the 
collision mechanisms and we briefly present the formal derivation of the SHE model. Then, 
we derive from this latter model the novel subband ET model. Finally, a subband drift- 
diffusion equation is obtained, as limit when the relaxation time goes to infinity. Section [3] is 
devoted to the numerical issues. Subsection 13. II presents the complete stationary model, the 
mixed finite elements scheme is described in Subsection 13. 2\ and the iterative approach is 
outlined in Subsection 13.31 Numerical simulations of an ultra-scaled Double-Gate MOSFET 
are presented in Subsection 13.41 

2 Formal derivation of adiabatic fluid- quantum models 
2.1 The quantum- kinetic framework 

We will assume in this work that the confinement direction is one dimensional whereas 
the transport takes place in a two dimensional domain. The domain is denoted Q = co x 
[0,£] with co C 1R 2 . The first two directions, called x G co, correspond to the classical 
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description of the gas, whereas in the third direction z G [0, £] quantum effects occur. The 
quantum confinement of the electron gas is described thanks to the eigen-elements of the 
ID Schrodinger operator. They are denoted (£ n , Xn)neN* and solve the eigenvalue problem : 

ft 2 d ( 1 d \ 

~\n I U ~X_ n ^-nXni 



2 dz\m*(z)dz' i fo 

Xn(x, •) G Hq(0,£), / XnXn> dz = 8 nn > . 

Jo 

In this equation fi is the reduced Planck constant, m* the effective mass. It is known that 
the eigenvalues € n in (12. ip form an increasing sequence tending to +oo. These functions 
depend on the potential energy defined by U = —eV, where e is the elementary charge and 
V denotes the self-consistent electrostatic potential, solution of the Poisson equation 

div x , z (e R (x,z)V x , z V) = -(N e - N D ). (2.2) 

Here Er(x, z) denotes the relative permittivity, Eq the permittivity constant in vacuum, 
Nd(x, z) is the prescribed doping density and N e (t, x, z) is the electron density. This density 
is described by a sequence of distribution functions {f n ) n &n* describing the repartition on 
each subband for the classical direction x G u> and the corresponding momentum variable 
k G M. 2 . It is written as 

+ 00 y „ 

f n (t,x,k)dk ) \xn\ 2 (t,x,z). 



N e (t,X,z) = / fn(t,X,k)dk\ 

n=l ' 



The evolution of distribution functions is governed by classical transport model in the x 
direction parallel to the gas. The total energy of the nth subband is defined by 

e n (t,x,k) = ^p^ + € n (t,x). (2.3) 
2m* 

Therefore the energy-band diagram of the semiconductor crystal is spherically symmetric 
and strictly monotone with respect to \k\. Then the Brillouin zone (which represents the 
elementary cell of the dual lattice L*) is equal to M 2 . Moreover, we point out that, in 
contrast with the classical counterpart, the energy-band depends on space and time. In a 
kinetic collisional framework, the distribution function f n of the nth subband satisfies the 
rescaled Boltzmann transport equation jU [31] : 

2 

a 2 d t f: + a{V k e n ■ V,/ n Q - V x e n ■ V fc / n Q ) = Qi d {f a ) n + jQe(f a ) n , (2.4) 

where Qzd is the collision operator for the lattice defect collisions and Q e is the collision 
operator for the elastic, non linear electron-electron collisions, accounting for intra-band 
scattering as well as for transitions between subbands. a and /3 are dimensionless parameters 
that satisfy a <C /3 <C 1. 

The main classes of lattice-defects that we shall consider are impurities and phonons [1] : 

Qld(f) — Qimp(f) + Qph(f)- 
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The elastic character of the impurity scattering leads to 

Qim P (f)n(k) := V / 3>j[^/(fc, k')S(e n (k) - E n >(k'))(f n >(t,x, k') - f n (t,x, k)) dk', 

n'eN* * 2 

where 5 is the Dirac measure and the dependence on t, x of e and <J>* mp has been omitted. 
The cross-section is assumed to be symmetric: <f> l ™£,(k, k') = &,™ P ri {k' , k). 
The electron-phonon collision operator is considered as 

QM)n(k)= E f KUk,k') 

([(iV ph + l)5(e n (A;) - e n ,{k') + a 2 e ph ) + N ph 5(e n (k) - e n ,{k') - a 2 e ph )}f n ,(k')(l - rjf n (k))— 
[(N ph + l)8{e n ,{k') - e n (k) + a 2 e ph ) + N ph 5{e n ,(k') - e n (k) - a 2 e ph )]f n (k)(l - v f n ,(k')))dk', 



where again, $^ n / (k, k') = $^ n (fc', fc), is the phonon energy, r] > is a dimensionless 
distribution function scaling factor and the terms < l—r)f n < 1 express the Pauli exclusion 
principle. N p h is the phonon occupation number, given by the Bose-Einstein statistics 

N ph = (> 2 Wfe^) - l) _1 , (2.5) 



(2.6) 



with Tl the lattice temperature and ks the Boltzmann constant. Formally expanding the 
phonon collision operator in power of a 2 , we get 

Q a ph {f) = Q ph M) + ^ 2 QXi(f) 

where Q^ h 1 is of order 1 when a goes to 0. 

The electron-electron collision operator is given by j33] 

Qe(f)n(k) = 

E / K<r,s( k > k 'i k 'l) 5 ^ + £ nM - < - + h - k' - k[) 

[frf'sA 1 ~ lf»)Q ~ iM " " ^X 1 " ^l)] dk'dkldk[. 

The notation /n',i, and stands for f n '{ki), f r {k') and fs(k[), respectively. 
We define then the elastic collision operator 

Qo(f)n — Qimp(f)n + Qph,o{f)n 

= E / WM*) - e nl (k'Wn>(k') - /„(*)) dk' , 

where $° , = « + (2iV p , + l)< n ,. We set 



(2.7) 



n,n' 

1 



QUf) = QphAf) + -Mf)- (2.« 



Then, the kinetic equation, starting point for the diffusive limits, is written in the fol- 
lowing form 

a 2 d t f% + a(y k e n ■ V x f* - V x e n ■ V fc /„ a ) = Q (f a )n + a 2 Q«(f a )n- (2.9) 
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2.2 Definitions and notations 



We first recall the coarea formula : for any C 1 function g : B (—>■ R, and any test function 
ip G C°(B), we have : 

where dS e (k) denotes the Euclidian surface element on the manifold g~ 1 (e). We denote 
dN e (k) = dS £ (k) /\S7 g{k)\. Taking g{k) = \k\ 2 h 2 /(2m*), the set of possible wave vectors of 
electrons belonging to the n — th subband and having total energy e is given by S £ ^ tn = 
{k G M 2 s. t. \k\ 2 = 1m*fr 2 {e - € n )} and dN £ _ tn = rf§2^%; where dS e - £n is the surface 
measure of the ball S £ - f . The coarea formula leads to : 



n r+oo / 

/ Mk)dk= / 



neN* 

and 

|2fe2 



ip n (k)dN e _ en (k) de, (2.10) 



£ / Mk)S(€ n + l -^-e)dk=J2[ ^ 



(k)dN e . en (k). 



ngN* K ngF 

Definition 2.1 We wn'W use t/ie following notations : 

• The density of states is defined by : 

N(t,x,e) := f dN £ . tn (k) = 2n^Af(t, x, e), 

where Af(t,x,e) = max{n GN* / 6 n (t, cc) < e} ; imi/i £/ie convention of Af(t,x,e) = 
i/e < €i(t,x). 

• TTie Fermi-Dirac function is given by 

^(t, x, e) = (rj + exp (^^r^ 

• We introduce the two Hilbert spaces 

+oo „ 

L 2 :={/ = (/nW, YW |/ n (fc)| 2 *<+oo}, 
endowed with the natural scalar product 

(f,g) f n (k)g n (k)dk, 

and 

L 2 . = {/ G L 2 (M) s. t. [ f 2 (e) d£ < +oo}, 



endowed with the weighted scalar product defined by 

de 



(f,9)r= / f(e)g{e 



T(e)(l - V T(e)Y 
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We will make the following assumption on the cross-section : 
Assumption 2.2 The coefficient $° n , satisfies, for Xq and X± two positive constants,' 

< A < K,n>N(t, X, En) < Al, $° n ,n>(k, k') = K>,n(k', fc), 

where N is the density of state defined above. 

2.3 First macroscopic scaling : the spherical harmonic expansion 
model 

For the sake of completness of this work, we present in this section the limit a — > of the 
kinetic equation ( 12.91) . All calculations will be done formally and we refer the reader to [7] 
where the rigorous derivation is studied. We consider the Hilbert expansion 

f a = f° + af 1 + a 2 f 2 + --- 

By linearity of the operator Q and by identifying the term of equal powers of a in ( 12. 9ft . 
we obtain 

Qo(/°) = 0, (2.11) 

QoifX = ^ • vj n ° - v,e n • v k fl (2.12) 

m* 

QoU'X = d t f° n + ^ • V x f l n - V x € n ■ V k f l n - Ql(f%, (2.13) 

m* 

where Q\ is obtained by taking a = in the expression (12.81) . 

We will then make use of the following properties of the collision operator (see [7]). 

Proposition 2.3 Under Assumption \2.2\ the elastic collision operator Qq defined in ( 12. 7h 
satisfies the following properties : 

1. The linear operator Qq : L 2 y L 2 is a bounded, symmetric, non-positive operator. 

2. For any bounded function ip : JBL \— > K., we denote ip(e) n {k) = ip( ^-^ + 6 n ). Then, 

V/a 2 , Qo(^(e)f) = ^(e)Q (f). 

3. The Kernel of Qq is the set 

Ker Q = {f e L 2 , s.t. 3ip : R -> R, / = ^i(e)}. 

^. T/ie range R(Qo) is closed and coincide with the orthogonal of the kernel of Qq given 
by : 

(Ker Qo) 1 - = {f e L 2 , s.t. ^ / f n (k) dN £ . 6n (k) = 0, for a.e. e > d}. 
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From Proposition 12.31 and (12.11)1 . we deduce that f° is an energy dependent function : 

f°(t,x,k) = F(t,x,e n ). 

Choosing ip : R — > R such that kip(e) G L 2 , we deduce from Proposition I2.3l that there exists 
a unique solution £ in (Ker Qq) 1 - such that 

-Qo(0 = ^(e). 
m* 



We can write £ = • ip(e) and from the second item of Proposition 12. 3[ is independent 
of the choice of the function ip. Then, equation (12. 12[) leads to 

ti(t, x, k) = -e n (t, x, k) ■ V x F(t, x, e n ). (2.14) 

Finally, the solvability condition of equation (I2.13P is that the right hand side belongs to 
(Ker Qo) ± - This leads to 

Yl I Wn + ^kS n ■ V x f l n - V x € n ■ V k f n - QlUX) dN £ _ en (k) = 0, for a.e. e > 6 1 . 

neN* S £-£n 

(2.15) 

Let us denote 

S e (F)=J2 f Q e (F) n dN £ _ en (k), and S 1 (F) = ^ / Q° ph>1 (F) n dN £ _ en (k). 

(2.16) 

Multiplying (12.151) by an energy-dependent test function 0(e) and integrating with respect 
to the variable e, we obtain for the first term : 



/oo n poo n 

/ &tfn dN £ _ 6n (k)(f)(e) de = J2 (d t F + d £ Fd t 6 n ) dN £ „ £n (k)(j>(e) de 

/oo poo I r 

Nd t F<f){£) de+ d £ FlJ2 d ^ \ dN e-eM \ 4>{s) de. 



Using the coarea formula (I2.10p . we deduce that 



/(E / (V k e n -VJ^-V x e n -V k ti)dN^ en (k)U(e)de = 

VneN* JSz-tn J 

= E / (v.-(|?^)-Vj fc .(/ 1 iv se e B M eB )dfc = v,.[5] / ^me n )a 

neN* jR2 m VneN* 7lR2 ™ 

f 00 r h 2 k \ 

= ~ V * [J2 ^ ® & n dN £ „ en (k) ■ V X F <f)(e) de, 

where the last identity is a consequence of (12.14ft . We define the diffusion matrix by 

D(t,x,e):= V / ^^e n dN £ _ en (k), (2-17) 
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and the current density by 

J(t, x, e) = -D{t, x, e) ■ V x F{t, x, e). (2.18) 

With these notations, we get that in the distributional sense, equation ( 12. 15ft is equivalent 
to the spherical harmonic expansion (SHE) model 

Nd t F + V x ■ J - Kd £ F = ^S e (F) + S\(F), (2.19) 



where k is given by 



x, e) = -2?r^ ( ( £ - e n) + ) • (2-20) 



The notation u + = max{0, u} denotes the positive part of u. We recall moreover a property 
of the diffusion matrix D stated in Lemma 2.8 of [TJ. We point out that the effect of the 
confinement is reflected in the special form of the coefficients of ( 12 . 1 9j) . which involve the 
subband energies. 

Lemma 2.4 The diffusion matrix D(t, x, e) define in ( 12. 17ft is a symmetric and nonnegative 
2x2 matrix. 

Remark 2.5 We end this section with a particular choice of the cross-section $° n ,, which 
allows to compute explicitly the diffusion matrix. If the cross-section is an energy- dependent 
function of the form 

<f>° ntn ,{t,x,k,k') = $°(t,x,e n ), 
then, after a straightforward computation, we have that 

Q n (t,x,k) = . . . - Vfc£ ra - 
<f>°{t,x,e n )N(t,x,£ n ) 

Therefore, the diffusion matrix defined in ( I2.17P has the expression 

TTh 2 /m* sr^ f \ -j- .X 



nGN* 



2.4 Second macroscopic scaling : the energy-transport model 

We start from the SHE model (12.191) and we assume that the electron-electron collision 
operator is dominant with respect to the second order correction of the phonon collision 
operator and therefore <C 1, in order to obtain an ET model. Passing through the 
SHE model, instead of starting directly from the Boltzmann equation allows to get an 
explicit expression of the coefficients, which is needed for numerical purpose. Moreover, 
the considered dominant scattering mechanisms provide an energy relaxation term in the 
macroscopic limiting model. 

The formal limit /3 — > in (12 . 19[) is again performed by means of a Hilbert expansion 

F = F° + /3F 1 + ■ — 
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Identifying equal powers of f3 implies 

S e (F ) = 0, (2.22) 

Nd t F° + V X J° - Kd £ F° - Si(F°) = DpoS^F 1 ), (2.23) 
where DpoS e denotes the Frechet derivative of S e at F°. 



We summarize below some useful properties of the collision operator S e defined in (12 . 1 6[) 
and of its Frechet derivative. 

Proposition 2.6 Under micro-reversibility assumptions on the cross-section <3> e , the oper- 
ator S e satisfies the following properties : 
(i) For all f,g G £ 2 (R), we have : 

[ S e (f)(e)g(e)de = 

~\ E / K.n',rAS k [f(Br(k'))f(e s (k[))(l - V f(e n (k)))(l - 
4 n,n>, r,s J 0P)* 

- f(e n (k))f(e n ,(h))(l - r]f(e r (k')))(l - vf(^(K)))[g' r + g' sA - g n - g n >,i] dkdk.dk 'dk[. 
(ii) kernel : 

Ker S e = {f e L 2 (R)-3fi(t,x),T(t,x) s. t. f{t,x,e)=^ >T {t,x,e)}, 
where J~h,t is the so-called Fermi-Dirac distribution function (see Definition \2.1\) . 

Proposition 2.7 The linear operator DjrS e satisfies 

(i) DjrS e is bounded, symmetric, non-positive on L^r. 

(ii) The kernel of DjrS e is given by 

Ker (D T S e ) = Span{F(l - rjF), F{1 - 7]F)e} 
(Hi) The range of DjrS e is closed and we have 

R(D T S e ) = Ker (D T S e ) L = {/ e L% / jf f(e) ( * ) de = J . 

These properties are an easy consequence of Proposition 3.16, Proposition 3.17 and 
Proposition 3.19 of [7j, using the fact that : 

r r+oo r 

/ S e (f)(e)g(e)de = £ / Q e {f) n {e)g{e) dN e _ en (k)de 
= H / Qe(fUe)g(e)dN e _ en (k)de 



/ Qe(f)n(£n)g(£n)dk, 
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where we use the coarea formula ()2.1(jp for the last identity. 

Formal derivation of energy-transport model. Let us come back to the formal limit 
(3 — > in f)2.19p . Thanks to Proposition I2.6[ equation ( 12 .22 j) implies that there exist p(t,x) 
and T(t, x) such that 



F°(t,x,e) = 7^ T (t,z,e) 



(2.24) 



From Proposition 12.71 we deduce that equation ( I2.23P admits a solution iff 

J (Nd t F° + V X -J°- nd £ F° - S^F )) ( M de = 0. (2.25) 

For the first term, the definition of the density of states N (see Definition 12. ip implies that 
N(t,x,e) = 27r^2-n if e 6 [£ n , € n +i) and vanishes for e < 6\. Then 



Nd f F° 



de 



E 



Tfl 

2ir—nd t F 



de 



Using the expression of the current ( I2.18p . we can rewrite the second term of ( 12. 25ft : 
From ( I2.20p . we deduce 

= - E 2^a,e„F»(e„) ( ^ ) - E 2l 4^» />H I J 



where we use an integration by part for the last identity. Finally, the solvability condition 
( I2.25P writes in the following form : 



0, 
-V* • 



NF 



de 







J^D(t,x,e)-V x F°(t,x,e) (^de = J^F ) (^de. 



(2.26) 



Let us denote by p and p£ the charge density and the energy density, respectively, 
associated to the Fermi-Dirac distribution function J-^t '■ 

p^ T (t,x):= / NJ r ^ T (t,x,e)de = ^2 / ^,T(t,x,e n ) dk, (2.27) 
Jr neN , </r 2 

p£p )T {t,x) := / NJ r ^ T (t,x,e)ede = / £ n J>,t(^ z, e n ) dfc. (2.28) 

^ „gN* </lR 2 
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We can state easily that for a Fermi-Dirac function, we have 

V x J>,r(*,a;,e) = -^>,t(1 - V^t) (eV x (-z—=) - V x ( 



/c B T k B T' 



Then, equation (I2.26P reads, using J R Si (J 7 ) cfe = 0, 

dtP^T ~ V x . • J x = 0, (2.29) 

3s (pS^t) + J2 27r Jp d t 6 " J ^(t, x, e) de-V x -J 2 = W, (2.30) 



where we denote 



J i: = / D(t,x,e)-V x ^ T (t,x,e)de = B 00 V x (-^)-B 01 V x (-^-), (2.31) 

JR KB J- KB-L 

J 2 := I eD{t ) x,e)-V x ^, T (t,x ) e)de = B 10 W x { 1 ^)-B 11 V x { T ^—) i (2.32) 

and where the diffusion coefficients are defined by 

%(t,a;) = / J D(t,x,e)£ i+ ^ iT (l - tte, for z, j = 0, 1, (2.33) 

with .D being defined in (I2.17p . The system (I2.29p - ()2.33p forms the energy-transport model 
in the transport direction for a partially confined electron gas. We recover the general form 
for energy-transport models for semiconductors (see |3j HI El] and references therein). The 
right hand side of the energy equation ( I2.30p is the so-called relaxation term W and it is 
defined by 

W:= f eS^^de = £ f + °° ( f Q° phA (f) n dN £ _ en (k)\ e de. (2.34) 

As for the SHE model (I2.19p . the presence of the subband energies in the diffusion matrix 
and in the relaxation term reflects the effect of the confinement in the transport equation. 

Then we have the following important property for the diffusion matrix which is an easy 
consequence of expression (I2.33P and Lemma 12.41 

Lemma 2.8 The diffusion matrix 

V := 



^00 ^01 

Kio On 



is a symmetric and positive definite matrix. 

Remark 2.9 After straightforward calculations, we can have an explicit expression of p and 
pS . In fact, 

. . ^ f En+1 m* ^_ 2nm*k B T ^ ,/i-€ n ,. 
Pn,T\t,x) = / 2n—nJ : fM:T (t,x,E)de= — ^ log(l + r/exp(— -— )), 
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p£^ T {t,x) = ^2 2n ~fi2 n / eJ c fJltT (t,x,e)de 

ngN* 



£n + l 



if rj > 0. For Boltzmann statistics rj = 0, we have 



P ,, T (t,x) = £ 2 7 r^T(M)exp( ^^ r( ^|' x) ) ) (2.35) 



ngN* 



p£^(t 7 x) = ^%k B T(t, x)(6 n (t, x) + T(t, x)) exp( //( ^ X) e " (t ' x) ). (2.36) 
2.5 Relaxation term W 

The relaxation term denned in (12.341) measures the influence of the interaction of phonons 
with the charge carriers. In [3]] we have formally established that a diffusion limit of the 
kinetic Boltzmann transport equation coupled to subband model in the scaling of domi- 
nant phonon-electron interaction leads to a drift-diffusion system in the transport direction 
coupled to the subband model. Using the coarea formula we can rewrite (I2.34p as 



Moreover, we have 



W= I e n QXi(^,T)dk. 



£ / e n Q a ph (^, T ) n dk = J2 [ AO 

nt^* J ^ 2 n,n> J ^ 

([(Np h + l)5(e n - e' n , + a 2 e ph ) + N ph 5(e n - e' n , - a 2 e ph )}e n J : '! n ,(l - rjTn)- 
[(N ph + l)5{e' n , - e n + a 2 e ph ) + N ph 5(e' n , -e n - a%/ l )]e n J r n (l - r)F' n ,)) dkdk', 

where the notation e' n ,, T n and T' n , stands, respectively, for £ n /(k'), J-^T^n) and J-*^^,). 
Since the Fermi-Dirac distribution function is energy-dependent, we have 

[ e n Q% h (^ T ) n dk = J2 [ $^,(A;,A;')(l-^)(l-^n) 

ngN* n,n' ^ 

([(JVph + l)«f(e n - + a% h ) + N ph 5(e n - e' n , - a 2 e ph )]e n M n ,- 
[(N ph + l)S(e' n , - e n + a 2 e ph ) + N ph 8{e' n , - e n - a 2 e ph )}e n M n ) dkdk', 

where M n = e^- £ ^)/( k BT) j g ^ e Maxwellian and it satisfies T n = M n (l — r/J-'n). Moreover, 

[ SnQ a ph {F»,T)ndk = Y J I KU^m-vK^-VJ'n) 

ngN* jR2 n,n> J ^ 

(5(e n - e' n , + a 2 e ph )[{N ph + l)M' n , - N ph M n ]e n - 
6(e' n , -e n + a 2 e ph )[(N ph + l)M n - N ph M n ,}(e' n , + a 2 e ph )) dkdk 1 . 
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By changing the variable e n with e' n , in the first term of the sum, we notice that only the 
a 2 -factor term does not vanish in this last identity, which can be rewritten as 

a 2 e ph )[N ph M' n , - {N ph + l)M n ] dkdk' = 
a% h )N ph M n (e »n r - e^ h ) dkdk', 

n,n> UM ' 

where the phonon occupation number N p h is defined in (12. 5p . Letting a — > 0, we have that 
N ph {e a2£ ^l kBT - e Q2 W fc ^) T L (i - ±). Thus 

W = T L d ~^r)J2 [ <n>(k, k')JF n {\ - r,Fn)e ph 6(e n - e' n> ) dkdk'. (2.37) 

The following lemma proves that If is a temperature relaxation term which relaxes T 
to the lattice temperature Tl. 

Lemma 2.10 Let W be defined in (I2.34p . Then, we have 

W-(T-T L )< 0. 

Proof. The proof of this result is an immediate consequence of (I2.3T[) . □ 



E / O 1 - - V^n)e ph S(s n - e' n , - 

n,n> 4 

E / <n'(* ~ VK>)(1 ~ r)T n )e ph 5{e n - e' n , - 



2.6 Formal derivation of drift-diffusion equation. 

In the case where the electron-phonon scattering is dominant, after a rescaling we have 
W = with a parameter 7 < 1, Then equation (I2.30P with expression (I2.37P implies 
that formally in the limit 7 — > 0, we have T = Tl- Then (I2.29P leads to the well-known 
drift- diffusion model 

d t p, - V, (l^ V ^) = °- ( 2 - 38 ) 

Moreover, assuming rj = 0, the equilibrium is then given by a Boltzmann statistics. We 
deduce therefore from (I2.35P that for T = T^ 

P, = 2^k B T L e^ TL E e~^ M . (2.39) 

nGN* 

Then 

We can introduce as in [8] the effective potential energy defined by 

V, = -k B T L \og ( E e~ enlkBTL ). (2.40) 

ngN* 
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Then, we have 

P/xV^/i = {k B T L V x p ll + Pfj,V x V s ) , 

such that, denoting D = 3 00 /p^, we recover from (I2.38P the standard formulation of the 
drift- diffusion system [8] : 

d tPtl - V x ■ (3(k B T L V xPfl + Pl N x V s )) = 0. (2.41) 

Remark 2.11 We conclude the section noticing the similarities and the differences between 
the classical and the partially confined energy-transport model. Due to partial confinement 
in the z direction, the electron density in the transport direction contains T as factor rather 
than T 3 / 2 . Moreover, we point out that the system (I2.29p - (l2.33p in the variables p/k^T, 
— 1/kjfT is in symmetric form, with the electric forces appearing in the diffusion coeffi- 
cients through the eigenenergies 6 n (which in turns depend on the electrostatic potential). 
In the partially confined framework, the diffusive limiting process brought directly to a set 
of variables which can be interpreted as dual entropy variables (using the denomination of 
non equilibrium thermodynamics fT% \Wj)- In this framework, T can be interpreted as elec- 
tron temperature, and the variable \i as quasi-Fermi potential energy. This fact is clear 
considering the drift- diffusion model (I2.4ip . Indeed, using the effective potential energy V s 
defined in f)2.40p . we can write the electron density (I2.39P as 

p, = 2^k B T L e^' k - T \ 

which gives the classical relation between electron density and quasi-Fermi energy. A chem- 
ical potential can then be defined as p c hem — A* — 

3 Numerical simulation 

The device we are modelling in this work is a nanoscale Double-Gate MOSFET (Metal 
Oxide Semiconductor Field Effect Transistor) such as in [31]. This device consists of a 
silicon film, characterized by two highly doped regions near the Ohmic contacts (denoted 
by source and drain) and an active region, called channel, with lower doping. The silicon 
film is sandwiched between two thin layers of silicon dioxide Si02, each of them with a gate 
contact. 

We assume invariance in the X2 direction (infinite boundary conditions), so that the 
problem is studied in a (xi, z)-domam. The device occupies a region of a 2D domain 
denoted by f2=[0, L\ x [0, £]. A schematic representation of the device is shown in Figure [TJ 

3.1 Energy-transport— Schrodinger— Poisson system 

In the following, we describe the collisional transport in the Double-Gate MOSFET, schema- 
tized in Figure [U by means of the energy-transport model (I2.29p -( |2.33p . The confinement 
is described by the subband decomposition approach, which involves the resolution of the 
eigenvalue problem (12. ip . taking also into account the presence of the oxide. Moreover, in 
order to provide explicitly computable diffusion and relaxation terms, the following physical 
assumptions will be used. 
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Figure 1: Schematic representation of the modeled device. 

Assumption 3.1 • (HI) The cross-sections $° and & ph are assumed to be energy de- 
pendent functions and to have the following expression (see /2l El [73]/ ) 

®l n ,{t,x,k,k') = <l?{t,x)e s n , $% n ,(t : x,k,k')=r%xK, (3.1) 

with —2 < s < 2. In t/ie physical literature, the values s = and s = 1/2 nave &een 
used J7J |g} [|f . 



• (H2) Tne electron density and the energy are assumed to be given by non-degenerate 
Boltzmann statistics, i.e. t] = 0, as in (I2.35P . (I2.36P . 

Using Assumption (HI), we deduce that the diffusion matrix (I2.2ip has the following 
expression 

D ^ = r^m ^. {e ~ eXM ' (3 - 2) 

where J\f(e) = max{n G N* / € n < e} is the number of non-zero terms in the sum (see 
Definition 12. ip . Moreover, thanks to the coarea formula in ( I2.37p . we have 



rp r+oo 

W = -An 2 e ph (l - ^ 



which under Assumption (H2), reads 



r+oo 

/ (f) ph e s Af 2 (e)J r (e)(l - ^(e)) de, 

J El 



+ 00 



W = -47T 2 e ph (l - ?±)e^ k ^ I <ff h e s Af 2 (e)e- £ /^ de. 

J ei 

By defining 

r+oo 

W = 47r 2 e ph / (j) ph e s M 2 {e)e~ e/{kBT) de, (3.3) 



we have the compact expression 



W = -W [l-^)e^ kBT K (3.4) 
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Assumption (H2) implies also that the density of charge carriers N e (t,x,z) is given by 
f -r~ , \k\ 2 h 2 \ 7 ,9 2irm*kBT v-^ <„ f \n- n 19 

N e = E / ^ft*. ^ + ^) rf/c|Xn| 2 = E e ( ^" )/fcflT |Xn| 2 . 

ngN* ^ M2 ngN* 

Finally, the coupled subband energy-transport model under Assumption 13. II is given by : 
Find fi(t,x), T(t,x), (6 n (t,x),x n (t,x)) for n > 1, and 2) such that 

dtP^-V x - J x = 0, in(0,L) (3.5) 
^(p^,t) + E ^^d t 6 n k B T(t, x)e^-^' kBT -V X -J 2 = W, in (0, L) (3.6) 

ngN* 



X n (t,a;, ■) G ^0(0,^), / XnXn'dz = 5 nn > 

Jo 



(3.7) 



div«,,( ea V^V) = - ( 2 -^L £ e^ lk ° T \Xn\ 2 -N D ), in fl, (3.8) 

where the expressions of and E^t with respect to the unknows are given in (12.351) - 
( I2.36p . In (13. 7h the effective mass m* takes different values in the Si and in the S1O2 
domain. Moreover, V c represents a given potential barrier between the silicon and the 
oxide. The currents Ji and J2 are given by the expressions 

Jx = ©ooV x .(/-) - D 01 V X .(-^-), (3.9) 
x k B T J k B T 

J2 = BioV sb (t^=) -BuV^r^-), (3.10) 
&Bi Kb-' 



where, under Assumption 13.11 the diffusion coefficients are given by 

AT{e 



p. . = J_ v / + °° ^ £ ^ e u^e)i*BT d 

*0 / A/" f c\ 6 



n6N* tn 

The relaxation term W is given by H3.4j) . 

This system is complemented with initial and boundary conditions. In particular, at the 
ohmic contacts and at the gate (see Figured]), we will impose Dirichlet boundary conditions 
for the potential, otherwise we fix homogeneous Neumann boundary conditions, which model 
isolating conditions. 

for z G {0,£},£ G Gate ; (3.11) 
for x G {0, L}, z G (0, £); (3.12) 

elsewhere, (3.13) 



V(x,z) 


— Vcate 


V(x,z) 


= v D , 


dV 
dv 


= 0, 
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where v is the outward unit normal. Since the transport occurs only in the longitudinal 
direction, we just have to impose boundary conditions in x = and x = L for \x and T . 
The temperature is assumed to be at the lattice temperature Tl, thus 

T(x)=T L , forx€{0,L}. (3.14) 

Then, we consider that the surface density of the charge carriers is almost constant near 
the frontiers x = and x = L and given by N%. The surface density being the integral over 
z of the total density (Njj = N + x we deduce 

/ N b h 2 \ 

fi{x) = fi b := k B T L log ' — _ e /k T , for x G {0, L}. (3.15) 

3.2 Stationary system 

Let us introduce the notations 

Then we can rewrite the expressions of the current fl3.9|) - fl3.10p as 

Ji^Boo^^V^ + Doi^^V^, (3.17) 
J2 = B w (u,v)V x u + B n {u,v)V x v, (3.18) 
where the diffusion coefficients are given by 

1 f +OQ e i+j ~ s (e — € ) 



hj{u,v) 



We define the relaxation coefficient in the same way : 



W (u,v) 



r+oo 

4ir 2 e ph ^ }h e s U 2 {e)e {e -^ v de. (3.20) 



Then, the stationary version of the energy-transport subband system ( I3.5l) -( l3~6l) in vari- 
able u and v writes 

-V x (B 00 (u, v)V x u + D 01 (u, v)V x v) = 0, (3.21) 
-V x (B 10 (u, v)V x u + Dn(u, v)V x v) = -W (u, v) (1 + k B T L v) , (3.22) 
Boundary conditions (13. 14[) — ( l3. 15[) become 

v(x) = v b := -7—=-, forxG{0,L}, 

and 

u(x) = u b : = log ' — — /k T , for x G {0, L}. 
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3.3 Numerical approach for the energy-transport system 

We introduce a partition of [0, L\ with nodes Xi, i = 0, • • • ,N X , and a partition of [0,£] 
with nodes Zj, j = 0, • • • ,N Z . We assume that the partitions are uniform and denote 
h = Xi — Then, we mesh the domain [0, L\ x [0, £] with rectangular triangles using the 
nodes (xi,Zj) previously defined. The Schrodinger equations and the Poisson equation are 
discretized with conforming P 1 finite elements. 

We consider here in details the discretization scheme for the equations governing u and 
v, assuming first that the eigenenergies € n are known. Using the following notations 

U = (u,v) T , J= (J 1; J 2 ) T , W(U) = (0 1 W (u,v)(l + k B T L v)) T 

equations (I3.17p - (l3.18p and ( I3.2ip -( )3.22p can be written in compact form as 

J = B(U)V X U, V X -J = W(U). (3.23) 

Denoting by Ui an approximation of U(xi), we take the piecewise constant approximation 
of U given, in the interval Ii := (xi-\,Xi), by 



2 

and define the piecewise constant diffusion coefficients and relaxation terms as 

B kl = B kl (U), for k, I = 0, 1; W= W(U). (3.24) 

We are going to use a mixed finite element discretization of lowest order in hybridized 
form [10] (see also [T9| [20] for applications to ET). Let us introduce the following finite 
dimensional spaces : 

X h = {a e L 2 (Q) : a(x) = a { + b^x - x Bi ) in I u i = 1, . . . , N x }, 
Yh = {£ £ L 2 (VL) : £ is constant in Jj, i — 1, . . . , iVx}, 

= {<? is defined at the nodes x Q , . . . ,x Nx , q(x ) = x(0), q(x Nx ) = 



where x Bi denotes the central point of the interval Ii, and \ is prescribed. 

Then, the mixed- hybrid formulation of f 1 3 . 2 3 j) reads as follows: Find G X%, G l^ 2 , 
and U h e Z h)Ub x Z/^ such that 

^^<f> h -B(U i y 1 J h dx+ ^Ph'Vsfadx-frh-Ud^ = 0, (3.25) 

Y^Q^ h -V x J h dx- Jw(JJi)-V h dxj = 0, (3.26) 



N x 



i=l 



18 



for all <frh G G Y^, and fi^ G Z% Q . Equation (I3.25P is derived from the weak 

formulation of the first equations in (13. 23ft ; ( I3.26P comes from the weak form of the second 
equations in (I3.23j) ; and finally, (I3.27P imposes the continuity of the currents at the nodes. 

Thanks to the discontinuity of the spaces Xh and Y"/,, we can apply static condensation 
in order to reduce the size of the discrete system and obtain an algebraic system for the 
variable Uh only. More precisely, choosing first the local basis 

(1,0) T in Ij f (0,1) T in J 4 

(0,0) T elsewhere ' 9h \ (0,0) T elsewhere 

in (I3.25P and then, analogously choosing ^ in (I3.26p . we obtain the piecewise linear (dis- 
crete) current 

J m = D(Z7 4 ) - + W(Ui)(x - x Bi ). (3.28) 

Imposing continuity at the nodes (through (I3.27P ) we obtain the final system 

- D(Fi)^_! + (B(Ui) + - B(U i+1 )U i+1 = -y (W(Ui) + W(U i+1 )), (3.29) 

for i = 1, . . . , N x — 1. We point out explicitly that, since the first component of W(U) is 
null, the approximation of the current J\ is piecewise constant (see ( 13. 28ft ) and that, thanks 
to (I3.27p . it is indeed globally constant. 

System (I3.29P forms a non-linear system in the unknown (u, v) that can be solved using a 
Newton algorithm. We point out that the Jacobian corresponding to this non-linear system 
can be easily computed noticing that, from the expressions ( 13.24ft and ( 13. 19p . we have 

dW M dW ki I-* dW M dW ke l-j 



du t dm-i 2 m dvi dv^ 2 fcm ' 
and that similar relations hold for the partial derivatives of W % . 

Remark 3.2 The expression (I3.19P is not practical for numerical purpose. However, with 
the definition N(e) = max{n G N* : € n < e} , we have 

™ n=lm=n Jern 



m=l u tm 



by interchanging the sums over m and n. We can rewrite 

1 

hj[u,v) 




e i+3-s e u+ev d 



m=l 



with the convention that Y^=i ^n/( m — 1) = for m — 1. Then, in the actual numerical 
computation, we can get an accurate approximation o/D^ by truncating the infinite sum to 
a finite number of eigenmodes. In fact, (6 n ) n forms an increasing sequence going to +oc, 
thus, since v < 0, we have that (e enV ) n fast decreases to 0. 
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3.4 Algorithm 



We are now ready to describe the algorithm used for the numerical resolution of the station- 
ary subband energy-transport model. The first step of the algorithm is the computation of 
the thermal equilibrium solution, with no applied drain-source bias. In this case the tem- 
perature and the Fermi level are constant along the device, therefore the problem reduces to 
solving the Poisson equation (13. 8p for a given temperature and Fermi level computed thanks 
to the boundary conditions. The computed potential at thermal equilibrium is used as a 
starting data for the following Gummel j21] iteration process : 

1. Lets V id be a given potential. 

2. We solve the eigenvalue problem (13. 7h on each vertical slice of the domain by diago- 
nalization of the Hamiltonian. Therefore we obtain the set {Xn{%i, Zj)} and {6 n (xj)}. 

3. We implement the Newton procedure which has been described above for the compu- 
tation of (u, v). 

4. We compute the density of charge carriers corresponding to the right hand side of 

(jSD 

Em* e u „. l9 
2vr— -e^Xn 2 - 
nr v 

nSN* 

We are then able to solve the Poisson equation ( 13.81) with boundary condition (13. lip . 
Indeed, due to the strong coupling of the entire system, the simple resolution of equa- 
tion (13. 8p does not provide a converging algorithm. Following [12J a Gummel iteration 
algorithm is used, amounting to compute the new potential V new by solving the fol- 
lowing modified Poisson equation 

dw(e R VV new ) + -N e ^p = - ( N e (l - ^) - N D ) , in SI, (3.30) 

"ref e V "ref J 

with V ref = k B T L /e. 

5. We repeat the last three steps until the quantity ||V^<z — Kieu)|U°° becomes sufficiently 
small. Once the convergence is reached, we increment the applied drain-source bias 
Vds °f 0.02 V and start a new iteration. 



3.5 Numerical results 

In this section we present and comment the performed numerical results. The modeled 
device is schematized in Figure [1] and the physical values are chosen as the ones in |3l] 
which are recalled in Table [TJ We take N x = 50 points in the transport direction and 
N z = 50 in the confined direction for all the tests. The results are presented for s = \ in 
(I3.ip . which corresponds to the so-called Chen model |14) . 

Some other physical coefficients should be determined. The effective mass is m* = 
0.19 m e with m e the electron mass, the lattice temperature is = 300 K and the scattering 
coefficient for the elastic collisions is the one used in |31| 0° = -5^7, where the low field 
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Table 1: Table of the main physical values 



Parameter 


Value 


Length 


Value 


N + 


10 26 m -3 


L s 


lOnm 


N~ 


10 21 m" 3 


L c 


30nm 


u c 


3 eV 


L D 


lOnm 


e R [Si] 


11.7 


P 

""OX 


3nm 


e R [S l 2 ] 


3.9 


isi 


5nm 



mobility is taken as fj, = 0.12 m 2 V~ 1 s~ 1 and the intrinsic density is given by n» = 10 10 m~ 2 . 
We have to fix the value of the scattering coefficient (p ph for the electron-phonon interaction. 
As noticed in subsection 12. 6[ when <ff h — > +oo, the model converges formally to the subband 
drift-diffusion system presented in [8j [31]. It is then interesting to compare the numerical 
results for large and small values of <j) ph . Figure [2] displays the / — Vds characteristics 
for <p ph = lO _4 /0° and <p ph = 10 5 /</>°. As expected and as noticed in [2], the energy- 
transport model gives higher currents compared with the drift-diffusion model, due to the 
electron velocity overshoot within the channel. Figure [2] (right) shows the temperature for 
cff h = 1O 5 /0 O , confirming that we are in the drift- diffusion regime, where the temperature 
is constant. Indeed we check numerically that maxT — minT = 2.6584 10~ 7 K. 




Figure 2: Left : / — Vds characteristics for Vq = OV and Vq = 0.2V. The dashed line 
corresponds to <p ph = 1O _4 /0 O , the solid line corresponds to <p ph = 1O 5 /0 O which is a good 
approximation to the drift-diffusion model. Right : Temperature in the device for (f> ph = 
1O 5 /0 O . We see that the temperature is almost constant. 

In the rest of the section, we present the results only for cf) ph = 1O~ 4 /0 O which corre- 
sponds to the energy-transport regime. Figure [3] displays the computed current vs drain- 
source applied bias characteristics with this chosen value. We present in Figure [3] (left) the 
characteristics for different numerical values of the Gate voltage Vg and with l$i — 5 nm. 
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In Figure [3] (right), we display the characteristics for Vq — V and for different geometry of 
the devices : l$i = 4, 5 or 7 nm. These characteristics are comparable to the one obtained 
in [2j [301 |3T]. We present in Figure [4] the evolution of the temperature in the device with 
respect to the drain-source voltage for two different values of the Gate voltage. Figure [5] 
displays the evolution of the mean velocity, defined by J\j '(qp^r), where the one dimensional 
density is given in (12. 35ft and the current in (13. 17ft . As expected, we notice an overshoot of 
the velocity at the frontier between the channel and the drain for high value of Vds- 
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Figure 3: I — Vds characteristics for different Gate voltages Vq with £$i — 5nm (left) and 
for different width of the Silicon in the DG-MOSFET with V G = V (right). 
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Figure 4: Evolution of the temperature in the device for a Gate voltage Vq = 0V (left) and 
V G = 0.2V (right). 

We plot in Figures [6] the shape of the density in the device for two differents drain-source 
voltage. For Vds — 0V, we are at equilibrium and the density is symmetric in the device. 
For Vds — 0.5V, we notice transport of the charge carriers in the device. 
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Figure 5: Mean velocity for different drain-source potentials Vds an d for Vq = OV (left) 
and V G = 0.2V (right). 
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Figure 6: Density of electrons in the device for Vds — OV (left) and Vbs = 0.5 V (right); in 
this simulation we take Vq = 0.2V. 

4 Conclusion 

A coupled quantum-classical model has been obtained for describing the transport of a 
partially confined electron gas. In a subband decomposition framework, the transport model 
is obtained by means of diffusive approximation from adiabatic quantum-kinetic models. 
The final system in the transport direction is obtained through two steps. First, under the 
assumption of dominant elastic scattering, a SHE system is derived (referring to [7]). Then, 
under dominant electron-electron collisions, an energy transport model is given, obtaining 
diffusion coefficients well suited for numerical purposes and with a relaxation term taking 
into account the electron-phonon interactions. In particular, with energy dependent cross- 
section of the collision operator explicit expression of the diffusion coefficients and of the 
relaxation term is derived and used for the numerical simulation of transport in a Double- 
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Gate MOSFET. In the limit of large electrons-phonons collisions, we recover the model of 
|31) . We point out that a saturation of the current is observed without need of resorting at 
mobility modeling as done in |31| . 

Acknowledgements. The authors acknowledge partial support from the Galilee Project 
no 25992ND of the Hubert Curien program : "Modeles numeriques du transport collisionnel 
dans des dispositifs nano-electroniques". 
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